## use Ben's nature genetics paper SNP array data to validate BisSNP accuracy. It used to test BisSNP after each upgrade
## author: Yaping Liu  lyping1986@gmail.com

use strict;
use Getopt::Long;
use File::Basename;

my $omit_validation="";
my $omit_genotyping="";
GetOptions( "omit_genotyping" => \$omit_genotyping,
			"omit_validation" => \$omit_validation,
			);

my $jar=$ARGV[0] || die "need input BisSNP jar file\n";
my $result=$ARGV[1] || die "need output result fiel name\n";
my $confidence=$ARGV[2] || 20;

my $use_age = "USAGE: perl validateBisSNP.pl [options] BisSNP.jar result.log.txt [confidence_score]\n\n require 6 CPUs and 8G memory";
if($ARGV[0] eq ""){
	print "$use_age\n";
	exit(1);
}


#my $REF="/home/uec-00/yapingli/software/gatk_read/Sting/resources/Homo_sapiens_assembly18_without_random.sort.fasta";
my $REF="/export/uec-gs1/laird/users/yaping/data/genome_data/genome/hg18/hg18_unmasked.plusContam.fa";
my $DBSNP="/export/uec-gs1/laird/users/yaping/data/genome_data/dbSNP/dbsnp_135.hg18.sort.vcf";
my $TESTBAM="/export/uec-gs1/laird/users/yaping/code/test_bissnp/accuracy_data/normalMerge_1MSNP.nodups.calmd.5beofre.5after.sorted.header.bam";
my $LOC="/export/uec-gs1/laird/users/yaping/code/test_bissnp/accuracy_data/Normal_SNP.filterCNV.txt.minReadCov_1.snpcaller.bed";
my $SNPARRAY="/export/uec-gs1/laird/users/yaping/code/test_bissnp/accuracy_data/Normal_SNP.filterCNV.txt.minReadCov_1.snpcaller.txt";

my $cytosine_out="BisSNP_test.".basename($jar).".cytosine.vcf";
my $snp_out="BisSNP_test.".basename($jar).".snp.vcf";

my $cmd = "java -Xmx8g -jar $jar -T BisulfiteGenotyper -R $REF -I $TESTBAM -D $DBSNP -L $LOC ";
$cmd .= "-vfn1 $cytosine_out -vfn2 $snp_out -stand_call_conf $confidence -stand_emit_conf 0 -dt NONE  -bsRate 0.9975 -out_modes EMIT_VARIANT_AND_CYTOSINES -nt 6 ";
$cmd .= "-mmq 30 -mbq 0\n";
print STDERR "$cmd\n";
if($omit_genotyping eq ""){
	system($cmd) == 0 || die "not finished BisSNP test call SNP module!!\n";
}

if($omit_validation eq ""){
	open(OUT,">$result") or die "cant open result file: $result\n";
	print STDERR "Validate SNP calling accuracy...\n";
	&test_snp_accuracy($snp_out);
	print OUT "\n";
	print STDERR "Validate Cytosine calling accuracy...\n";
	&test_methy_accuracy($cytosine_out);
	close(OUT);
}


sub test_snp_accuracy{
	my $file=shift @_;
	
	open(FH,"<$file") or die "cant open BisSNP snp.vcf file: $SNPARRAY !\n";
	my @vcf=<FH>;
	chomp(@vcf);
	close(FH);
	
	open(FH,"<$SNPARRAY") or die "cant open validate SNP array data file: $SNPARRAY !\n";
	my @in=<FH>;
	chomp(@in);
	close(FH);
	my $count = 0;
		my $count_het = 0;
		my $count_array_het = 0;
		my $count_hom = 0;
		my $fp = 0;
		my $fn = 0;
		my $total_count=0;
		my $ag_count=0;
		my $count_array_hom = 0;
		my %vcf_hash = ();
		foreach my $line(@vcf){
			next if $line=~/^\#/;
			my @splitvcf=split "\t",$line;
			my $key = $splitvcf[0]."~".$splitvcf[1];
		#	next if $splitvcf[5] <=30;
			$vcf_hash{$key} = $line;
		}
		foreach my $line(@in){
			my @splitin=split "\t",$line;
			if($splitin[3] ne $splitin[4]){
				$count_array_het++;

			}
			else{
				$count_array_hom++;
			}
			my $key = "chr".$splitin[1]."~".$splitin[2];
			if( exists $vcf_hash{$key} ){
				$count++;
				my @splitvcf=split "\t",$vcf_hash{$key};
				if(($splitvcf[9]=~/0\/1/ or $splitvcf[9]=~/1\/2/ or $splitvcf[9]=~/0\/2/) and $splitvcf[6]=~/PASS/){
					$count_het++;
					if($splitin[3] eq $splitin[4]){
						$fp++;
					}
				}		
				else{

						if($splitin[3] ne $splitin[4]){
							$fn++;
							if($line=~/\d+\tA\tB\tA\tG/ or $line=~/\d+\tA\tB\tG\tA/ or $line=~/\d+\tA\tB\tT\tC/ or $line=~/\d+\tA\tB\tC\tT/){
								$ag_count++;
							}
							
						}
						$count_hom++;

				}
			}
			else{
				if($splitin[3] ne $splitin[4]){
					$fn++;
					if($line=~/\d+\tA\tB\tA\tG/ or $line=~/\d+\tA\tB\tG\tA/ or $line=~/\d+\tA\tB\tT\tC/ or $line=~/\d+\tA\tB\tC\tT/){
						$ag_count++;
					}
					
				}
			}
			$total_count++;
		}
		my $fp_rate = $fp/$count_array_hom;
		my $fn_rate = $fn/$count_array_het;
		my $err = $fp_rate + $fn_rate;
		my $fn_1_rate=1-$fn_rate;
		print OUT "category\terror_rate\tfalse_postive\tfalse_postive_rate\tsensitivity\tfalse_negative\tfalse_negative_rate\tarray_het_snp\tarray_hom_snp\thet_snp_by_bissnp\tct_snp_in_array\n";
		print OUT "SNP_accuracy\t$err\t$fp\t$fp_rate\t$fn_1_rate\t$fn\t$fn_rate\t$count_array_het\t$count_array_hom\t$count_het\t$ag_count\n";
	
}

sub test_methy_accuracy{
	my $file=shift @_;
	
	open(FH,"<$file") or die "cant open BisSNP snp.vcf file: $SNPARRAY !\n";
	my @vcf=<FH>;
	chomp(@vcf);
	close(FH);
	
	open(FH,"<$SNPARRAY") or die "cant open validate SNP array data file: $SNPARRAY !\n";
	my @in=<FH>;
	chomp(@in);
	close(FH);
	
		my $count_array_not_C = 0;
	
		my $fp = 0;
		my $fn = 0;
		my $count = 0;
		my $count_array_hom_C = 0;
		my $true_pos=0;
		my $true_neg=0;
		my %vcf_hash = ();
		foreach my $line(@vcf){
			next if $line=~/^\#/;
			my @splitvcf=split "\t",$line;
			my $key = $splitvcf[0]."~".$splitvcf[1];
			$vcf_hash{$key} = $line;
		}
		foreach my $line(@in){
			my @splitin=split "\t",$line;
			$count++;
			if(($splitin[3] eq $splitin[4]) and ($splitin[5] eq "C" or $splitin[5] eq "G")){
				$count_array_hom_C++;
			}
			else{
				$count_array_not_C++;
			}
			my $key = "chr".$splitin[1]."~".$splitin[2];
			if( exists $vcf_hash{$key} ){
				my @splitvcf=split "\t",$vcf_hash{$key};
				if(($splitvcf[6] eq "PASS") and (($splitvcf[9]=~/0\/0/ and ($splitvcf[3] eq "C" or $splitvcf[3] eq "G" )) or ($splitvcf[9]=~/1\/1/ and ($splitvcf[4] eq "C" or $splitvcf[4] =~ /\bC\,/ or $splitvcf[4] eq "G" or $splitvcf[4] =~ /\bG\,/)) or ($splitvcf[9]=~/2\/2/ and ($splitvcf[4] =~ /\,C/ or $splitvcf[4] =~ /\,G/)))){
					if(($splitin[3] ne $splitin[4]) or ($splitin[5] ne "C" and $splitin[5] ne "G")){
						$fp++;
					
					}
					else{
						$true_pos++;
					}
				}
				else{

						if(($splitin[3] eq $splitin[4]) and ($splitin[5] eq "C" or $splitin[5] eq "G")){
							$fn++;
						}
						else{
							$true_neg++;
						}
				}
			}
			else{
				if(($splitin[3] eq $splitin[4]) and ($splitin[5] eq "C" or $splitin[5] eq "G")){
					$fn++;
				}
				else{
					$true_neg++;
				}
			}

		}
		my $fp_rate = $fp/$count_array_not_C;
		my $fn_rate = $fn/$count_array_hom_C;
		my $err = $fp_rate + $fn_rate;
		my $fn_1_rate=1-$fn_rate;
		print OUT "category\terror_rate\tfalse_postive\tfalse_postive_rate\tsensitivity\tfalse_negative\tfalse_negative_rate\tarray_not_C\tarray_hom_C\ttrue_pos\ttrue_neg\n";
		print OUT "Cytosine_accuracy\t$err\t$fp\t$fp_rate\t$fn_1_rate\t$fn\t$fn_rate\t$count_array_not_C\t$count_array_hom_C\t$true_pos\t$true_neg\n";
}

